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Abstract 

A mixed basis approach based on density functional theory is extended to one-dimensional 
(ID) systems. The basis functions here are taken to be the localized B-splines for the two finite 
non-periodic dimensions and the plane waves for the third periodic direction. This approach 
will significantly reduce the number of the basis and therefore is computationally efficient for 
the diagonalization of the Kohn-Sham Hamiltonian. For ID systems, B-spline polynomials are 
particularly useful and efficient in two-dimensional spatial integrations involved in the calculations 
because of their absolute localization. Moreover, B-splines are not associated with atomic positions 
when the geometry structure is optimized, making the geometry optimization easy to implement. 
With such a basis set we can directly calculate the total energy of the isolated system instead of 
using the conventional supercell model with artificial vacuum regions among the replicas along 
the two non-periodic directions. The spurious Coulomb interaction between the charged defect 
and its repeated images by the supercell approach for charged systems can also be avoided. A 
rigorous formalism for the long-range Coulomb potential of both neutral and charged ID systems 
under the mixed basis scheme will be derived. To test the present method, we apply it to study 
the infinite carbon-dimmer chain, graphene nanoribbon, carbon nanotube and positively-charged 
carbon-dimmer chain. The resulting electronic structures are presented and discussed in details. 
PACS: 71.15.Mb, 73.20.-r 
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I. INTRODUCTION 


There has been increasing interest in one-dimensional (ID) systems on the nanoscale, 
such as tubes, wires, rods, ribbons, etc, because the electronic properties of these systems 
are fundamentally different from those in higher dimensions due to their unusual collec¬ 
tive excitations. Expectations concerning the creation and application of improved func¬ 
tional electrodevices with better performance characteristics are rising from the intensive 
exploration of ID systems. Particularly, the emergence of nanotechnology has led to the 
realization of ID materials and stimulated both academic research and material innovation. 

First-principles methods based on the density functional theory have proven to be pow¬ 
erful and successful in investigating the electronic structures and properties of solids. The 
use of a plane-wave basis is most natural for infinite 3-dimensional periodic systems, such 
as bulk solids, because of its easy implementation and the fact that the convergence of the 
calculation can be checked systematically. 

To retain all the advantages of plane-wave expansions and of periodic boundary condi¬ 
tions in the investigation of low-dimensional systems, which are hnite along the non-periodic 
direction(s), the conventional supercell approximation is adopted by introducing some artifi¬ 
cial vacuum space to separate the periodic replica along the non-periodic direction. However, 
this approach suffers from one main drawback. For example, in two-dimensional (2D) sys¬ 
tems, it requires a vacuum layer of large thickness such that the interactions between the 
adjacent slabs are negligible, and therefore increases the number of the plane waves along 
that direction. In particular, for charged systems (e.g. charged defects), the rather long- 
range tail of the Coulomb potential inevitably requires an extremely large separation of the 


two slabs and makes the calculation impractical [l|. Many correction schemes have been 
devised to remedy this difficulty Sj-jd]. These drawbacks become even more significant in 
ID systems. 

In previous work Sj-js], a mixed planar basis approach that is conceptually simple, has 
successfully been introduced for the first-principles calculations of 2D systems by expanding 
the wavefunction along the periodic directions with 2D plane waves but, for the finite non¬ 
periodic direction, with ID localized basis of Gaussian functions or B-splines {q]. The use 
of this mixed basis has several advantages over the supercell modeling: (1) It resumes the 
layer-like local geometry which appears in surfaces and describes the wavefunction in a 
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natural way. (2) Because one can calculate the total energy for an isolated slab instead of 
using a supercell consisting of alternating slab and vacuum regions, the physical quantity, 
such as the work function can be immediately obtained without any correction. (3) For 
charged systems, the spurious Coulomb interaction between the defect, its images and the 
compensating background charge in the supercell approach can be automatically avoided. 
(4) The number of the basis is signihcantly reduced, easing the computational burden for 
the diagonalization of the Kohn-Sham Hamiltonian. 

To preserve the above good properties, in the present work, we extend our earlier work 
8| to ID systems, i.e., with two sets of B-splines to expand the wavefunctions along the 
two hnite non-periodic directions and ID plane waves along the periodic one. B-splines are 
highly localized and piecewise polynomials within prescribed break points which consist of 
a sequence numbers called knot seq uence and have proven to be an excellent tool for the 
description of wavefunctions js]- 12|. B-splines have the following properties: (1) due to their 
absolute localization, the relevant matrices are sparse. This is particularly useful and efficient 
for ID systems, which involve multi-dimensional spatial integrations. (2) B-splines possess 
good flexibility to represent a rapidly varying wavefunction accurately with the knots being 
arbitrarily chosen to have an optimized basis. (3) B-splines are, unlike the atom-centered 
Gaussian basis, independent of atomic positions, so the geometry optimization can be easily 
implemented. Here, a rigorous formalism designed to treat the long-range Coulomb potential 
of both neutral and charged ID systems within the present mixed-basis framework is also 
developed. 

We apply the present approach to study the inhnite carbon-dimmer chain, graphene 
nanoribbon, carbon nanotube, and the case of positively-charged carbon-dimmer chains. 
We perform the band structure calculation using Vanderbilt’s ultra-soft pseudopotentials 
(USPP) [^. Extensive cornparisons are made to the standard supercell approach with 
the popular VASP code 1^, 15|. It is found that the calculated band structures are very 
promising but the number of the basis is signihcantly reduced. Aside from the reduction, 
no further corrections are needed for the charged chain. 
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II. METHOD OF CALCULATION 


A. B-splines 


For the sake of completeness, we first briefly summarize the B-spline formalism. More 


details can be found in Refs. 


^ and 


9|. 


B-spline of order k consists of positive polynomials of degree k — 1, over n adjacent 
intervals. These polynomials are determined by a knot sequence {rj} and vanish everywhere 
outside the sub intervals r* < s < Tj+k. The B-spline basis set is generated by the following 
relation : 


S-Ti 








with 


^*, 1 ( 5 ) = 


1, Ti <S < Tj+i 

0, otherwise . 
The first derivative of the B-spline is given by 


d 


B, 


K 




K 


■-Bi+l,re-l(s). 


( 1 ) 


( 2 ) 


( 3 ) 


ds Tj '^i+K hi+l 

Therefore, the derivative of B-splines of order k is simply a linear combination of B-splines 
of order k — 1, which is also a simple polynomial and is continuous across the knot sequence. 
Obviously, B-splines are flexible to accurately represent any localized function of s with a 
modest number of the basis by only increasing the density of the knot sequence where it 
varies rapidly js]. 


B. Relevant matrix elements within B-spline basis 

In Vanderbilt’s USPP scheme [^, the wavefunction 0* satisfies a secular equation of the 
form 

>= e,R|0, > (4) 

under a generalized orthonormality condition 

< (j)i\S\(j)j >= 6ij . (5) 

Here, 

H = -V^ + Vpp + Vh + W , (6) 
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with Vpp, Vhi and Vxc denoted as the pseudopotential, Hartree potential, and exchange- 
correlation potential, respectively, and S' is a Hermitian overlap operator. In the following, 
we will give a detailed description of the relevant calculations for H and S operators. 

By using two sets of B-splines to describe the non-periodic x and y directions and ID 
plane waves for the periodic direction, the present mixed basis used to expand (pi is dehned 
as 


< r| k-h G;j, k; j > = 


VI 




(7) 


where G and k denote respectively the reciprocal lattice vector, and the Bloch wave vector. 
L is the length of the system along the direction. 

We dehne 

< BaBe,.>.= jd,B, mb,A s) (8) 

where s = x or y. It is an integration of local polynomials with bounded support and 
vanishes unless the condition \i — i'\ < n is fulfilled . This property is particularly useful for 
higher-dimensional integrals involved in the calculations for ID systems. 

The overlap matrix elements between two basis states are given by 


< k + G]i, K]i , K I k -|- G';j, K;j , 


K 


> = < Bm\ Bi,, >x < 




Bi' ^1^1 


B 


j'M ^G,G' 


(9) 


The kinetic energy matrix elements are given by 

< k + G;i, K;i , K \ — | k -|- G'; j, K]j,K > (10) 

[ ^ B B j i^ ^ BaB jt f^f T ^ Bi^fi ^ B B j/p^f (71) 

+ < Bi^fi Bj^,^ >x < >y {k + Gfi ] 5 g,g' ■ (12) 

B'i^fis) is the derivative of Bi^,^{s). As mentioned above, because of the absolute local¬ 
ization of Bi^i^ and Bfi^, the evaluation for the kinetic part of H\pi > is only an order 

Nd{K-l)K, (13) 


where is the number of the basis. In practical applications, n is usually set to be 4 
or 5. Therefore, the computational effort for the construction of the kinetic energy matrix 
elements scales linearly as A^. 

The local part of Vpp on each atomic site with species a concerned here was htted as 


Id: 


loc 


r = - 


erf 


r 




(14) 


5 



Presently, the rf-like pseudopotential is chosen as the local pseudopotential. The hrst term 
on the right hand side of Eq. flTT|) will be referred to as the core term due to the core charge 
distribution 


nAr) = 


712 


_ 

e ^ 


The local pseudopotential of the crystal is then given by 


VLOc(r) = V;-(r - R") , 


(16) 


where R'^ denotes the position of each atom with species a. The Fourier transform of the 
local pseudopotential of the crystal excluding the core term, V-^qq, is given by 


^LOc(P) 9 ) 


-igz 



-alip-RZ? 


(16) 

(17) 


Here, is the length of the unit cell (UC) along 2 ; axis. 

The total charge distribution n is dehned as the sum of the core charge distributions for 
all atoms in the sample Uc, plus the electronic charge distributions ng. 


n(r) = nc(r) + ne(r). (18) 

The Coulomb potential due to the total electron charge distribution and the exchange- 
correlation potential 14c should be determined self-consistently. The exchange-correlation 


potential are deduced from the Monte Car 
and parametrized by Perdew and Zunger 17 


o results calculated by Ceperley and Alder [16 1 
. We write 


14c(r) ^ ^ Vxcjpi 9 ) 




(19) 


With the use of the ID Fourier transformation of 1/r (see the Appendix) 

- = - [ dq.KMp) 
r 71J 

where Kq is the modihed cylindrical Bessel function of zero order, the Coulomb potential 
due to the total charge distribution is given by 

n(r') 





dh' 


J2v^^\p,g) 


^igz 


( 20 ) 
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where 


with 


y^^\p,9) = ‘^ // d^p'nip',g)Ko{\g\\p- p'\), 


dz n(r) . 


( 21 ) 


n{p,g) = 


Kq will diverge when S' —)■ 0. But, 0) is still well dehned, as explained below. 


For the s 0 case, using the asymptotic behavior of Kq 

A''o(a;) —)■ — Inx + ln 2 — 7 , a; —)■ O’*" 

where 7 is the Euler-Mascheroni constant, 0) is now split into two terms, 

i/(^)(p, 0 ) = 2 (ln 2 - 7 -ln|c/|) [[ p'n{p',0) - 2 [[ d'^ p'n{p',0)\n{\p - p'\) . ( 22 ) 


The first term could be safely dropped if the system is charge neutral (ff d?p'n{p', 0) = 0). 
We demonstrate in the following that omitting such term still holds for charged systems. 

We assume L, the length of the system along the 2 ; direction, is arbitrarily large but hnite. 
With 

n{r) = ^ n{p,g) (23) 


the (7 = 0 component of will be 


vy\p,o) = 



LI2 


dz' 


n{p', 0 ) d'^p' 


(24) 


l-L/2 a/A p2 + ^'2 

where Ap = |p — p'|. Suppose L be much larger than the cell size in the xy plane, i.e., 
L Ap, then 

dz' , L /2 + VLV4 + Ap 2 


^/Ap' 


.2 _i_ 


= In 


Ap 


In L — In I Ap| 


So, 


Vy^(p,0) = 2lnL // d^p'n(p', 0 ) — 2 // d^p'n(p', 0 ) In 


IP-P 


(25) 


(26) 


In the case of charge neutrality, the first term on the right hand side of Eq. (l26ll vanishes and 
we retain Eq. fl 22 l) . For charged systems with net line charge density Ni = JJ d?p' n{p', 0 ) 7 ^ 
0, such term would be huge. But clearly it is a constant that is independent upon p, and 
only causes a shift to the total energy. This kind of constant is irrelevant to the band 
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structure calculation. Therefore, within the present mixed basis approach, we could, just 
like the charge-neutral case, omit the hrst term without further corrections for the charged 
systems. 

The numerical integration around the singularity aX p = p' in Eqs. fl^TD and ( [2^ can 
be carried out and averaged over one hner sub-grid unit. In practice, we perform such an 
integration over a circle C, whose area is equal to that of the grid unit, as shown in Fig. [1] 
We also assume n{p') be constant within such a circle. Then, 


2 11 cfp'n{p',g)Ko{\g\\p-p'\) ~ 47rn(p,^) 


4 - 

.9 \9\ 


(27) 


and 


d p'n{p', 0) ln(|p - p'l)-n(p, 0) 


c 


2 1 ^2 


( 28 ) 


Here, e is the radius of the circle C and Ki is the modihed Bessel function of order 1. 
We dehne the total local potential Ws as 

Wff = + w, + VI. 


XC w I'log • 


(29) 


To construct the >, we need to calculate < k-|-G; i, k; i', k |Kfr|0i > in the real- 

space. Again, it can be done efficiently with the advantage of the absolute localization of 
B-splines mentioned before. 

Now, lets turn to the nonlocal part of H. The atomic nonlocal potential in the Kleinman- 

n 

Bylander (KB) form [18| is given as 


Vff,{T-W) = 


n9r,R- >< 




nlmn'Vm' 


with 




•)0,(T 



(30) 


(31) 


The projector and the augmentation function Q^imn'i'm' vanish outside the atomic core 


region. denotes the projector centered at the R®" atom, i.e., 

Pndm 

The nonlocal pseudopotential of the crystal is 

rNL(r) = 5; - R') , 

cr,R‘^ 


o-.R"" 


R-"). 


(32) 
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In practical calculations, we define a box, which is large enough to contain the core region 
Q- The - Rpz) inside the box is transferred to G space, then 


KZ = - R") = E - «il • G) 


—iGRZ AGz 


G 


(33) 


Therefore, 


< k + G; i, k; i , K \(5, 


(7,R'^ 

ndm 


> 


= U‘pBUx)B,,y{y)P:jJp-RlG') dz . (34) 

JJ J 


Similarly, both Qnimn'Vm'i.^ Kff(r) in Eq- are also transferred by fast Fourier 

transfer (FFT), 


Q 




> - R") = 


Gh 


Q 




,(p — Rn^Gh] Rl) 


JGhZ 


(35) 


V,s{r) = J2Vesip,G^) 

Gh 


(36) 


Then, we obtain 




nMV = + LJ E // GOl* - -Rf. Gy, flj) 

Gh 


(37) 


Lj is the length of the core region box along the ^ axis. Note that the FFT grid density 
for Gh in the summation of Eqs. fl3^ and fl36|l is not necessarily the same with that for the 
wavefunction [19[. in the above should be calculated self-consistently. From Eqs. 


and (jST]), we can evaluate V^i,\cj)i >. In doing this, we need to calculate 


< k +G; j',^'IE nlI^* >= ^ ^ 


nlm^n'l'm' 


< k + G; j, k-j\k >< > 


TCr,R‘^ 


cr,R‘^ nlmn'l'm' 

(38) 

Because of the KB separation form in Eq. (I38|) . the computation effort for this part is also 
linear to Nd- 

As for the Hermitian overlap operator S', which is peculiar to the Vanderbilt USPP 
scheme, it is given by 


S 






(t,R'^ 

n,lm 


(TjU®" nlmn'l'm' 


X /5 


o-,R'^ I 

n' ,l'm' I ’ 


(39) 
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where q^imyi'm' = III S\(j)i > will be obtained similarly. Finally, the 

charge density from the wave function is augmented inside the core region, 


Peir) = +Y1 Qnlm,n'l'm'( 

i cr,R‘^ nlmn'l'm' 


R") < X > ]■ (40) 


cr,R“ 


To hnd the lowest eigenvectors of the Hamiltonian matrix H, we used the Lanczos-Krylov 
method developed previously js], with the orthonormality of the wavefunctions maintained 
throughout by the standard Gram-Schmidt orthogonalization procedure. 


III. APPLICATIONS OF PRESENT METHOD 


A. infinite carbon-dimmer chain 


For the hrst example, we study a simple testing case of inhnite carbon-dimmer chains, 
as shown in the top panel of Fig. |2j The chain has a periodicity of \/3ao (oq = 2.461 A) 
along the z axis and the atom distance of C-C is chosen to ao/\/3. Two sets of 11 B-splines 
that each is dehned over a range of 3.25 Oq, are used to expand the x and ^/-component 
wavefunction, respectively. The energy cutoff of the ID plane waves along the z axis is 20 
Ry. Special fc-points of 1/8 and 3/8 (in unit of 27r/\/3ao) were taken to sample the ID 
Brillouin zone. The C USPP was generated from the Vanderbilt’s code 2^ and its quality 
was examined previously Sj]. For comparison, we also performed the calculation by using 
the VASP code with the projector-augmented-wave potential (PAW) jl^. Q. A typical 
vacuum space of 10 A x 10 A required in VASP was used in the calculation. The potential 
is determined self-consistently until its change is less than 10“® Ry. 

Figure|2](a) displays the band structure between F (0) and X (1/2). The VASP counter¬ 
part is shown in Fig. [2](b). Clearly, the present calculation agrees nicely with that by VASP. 
The splitting of the twofold degenerate bands due to the symmetry of the system is found to 
be smaller than 0.0001 eV. Note that the number of the basis is significantly reduced from 
~ 4300 by VASP to ~ 1900 by the present method. 


B. graphene nanoribbon 

Next, in order to provide a stringent test of the present method, a realistic system of 
the armchair graphene nanoribbon is considered. The width of the ribbon here is chosen to 
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include Na = 19 carbon dimmer lines, as indicated in the left part of Fig. [3l One set of 
73 B-splines distributed over a range of 14.25 oq and another of 11 B-splines over a range 
of 3.25 Qq were used to expand the x- and ^/-components of the wavefunction, respectively. 


le periodic direction was 20 Ry. The atomic 
8|. The special fc-points for sampling the ID 


The energy cutoff of the ID plane waves for t 
positions of the system were taken from Ref. 

Brillouin zone are 1/8 and 3/8. 

The present band structures near the Fermi level are displayed in FD. |3l^a). To make a 
comparison, we also show the previous calculations and VASP results [Sj in Figs. Hl^b) and 
(c), respectively. In previous work js], we adopted the planar mixed basis set to study the 
nanoribbon by using the surface supercell modeling for the x — z plane. Namely, we used 
only one set of B-splines to expand the y component of the wavefunction. As compared 
with the VASP calculation, which was obtain by the traditional 3-dimensional supercell 
modeling, the total number of the basis in previous work js] is reduced from ~17000 to 
~ 12300. Now, if we use another set of B-splines for the x non-periodic direction, then the 
number is further reduced to 8800 only. However, it can be seen from the hgures that there 
is a very nice agreement between these three approaches. This reflects the advantage of 
using B-splines for non-periodic directions over the plane waves, especially for ID systems. 
Actually, we have done all the calculations of the Na < 19 families. All the results are found 
in excellent agreements with those by VASP. Therefore, we are convinced that the present 
program has been implemented successfully for ID systems and the results obtained are 
very reliable. Moreover, as compared to the traditional supercell modeling, the use of the 
present mixed basis will signihcantly reduce the number of the basis functions and speed up 
the calculations for ID systems. 


C. zigzag carbon nanotube 

Now, we study carbon nanotubes (CNTs), allotropes of carbon with a cylindrical nanos¬ 
tructure. We choose the (4,0) zigzag CNT (Fig. 0]) as the third example. 

A mixed basis set with two identical sets of 29 B-splines over a range of 6.25 Oq along the 
two non-periodic directions and the plane wave cutoff of 20 Ry along the periodic direction 
are used. All C atoms were kept at the ideal positions which were obtained by rolling the 
ideal graphene sheet with the C-C bond length taken to be ao/\/3. The special fc-point for 
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sampling the ID Brillouin zone was 1/4. The resultant band structures are presented in Fig. 
m^a), along with the VASP calculations in Fig. IH^b) for comparison. 

Clearly, all bands compare favorably with those by VASP. It is worth mentioning that it 
would be more efficient to investigate nanotubes or rods if we expand the wavefunction in 
cylindrical coordinates rather than in Cartesian ones. Here, we utilize this system to test 
the present program. With the present algorithm to study CNT, the B-splines used are 
more densely distributed as compared to the hrst two cases. And the relevant real-space 
integrations were carried out with hner grids (48 equal divisions within oq) to obtain precise 
results. Nevertheless, the total number of the basis used here is reduced from ~9700 by 
VASP to ~8000 by the present mixed-basis approach. To sum up, we have demonstrated 
that the present program is computationally efficient to produced reliable results. 


D. charged carbon-dimmer chain 


Finally, we apply the present method to charged systems which are very challenging 
for the supercell modeling because of the spurious long-range Coulomb interaction between 
the defect and its periodic images. For simplicity, we still use the same system of the 
inhnite carbon-dimmer chain, but here with one of every eight electrons removed. That 
is, the nominal ionicity of C in this artihcial positively-charged chain is -1-0.5. All other 
computational conditions are similar to the hrst example described above. The resulting 
band structure is shown in Fig. [5](a). We also show the VASP result in Fig. |5](b), which 
was obtained by using a homogeneous compensating background charge. 

It is obvious that the present approach yields very similar results to those by VASP. 
Notably, the convergence rate of the calculation is fast and comparable to the neutral case. 
In the supercell approach, the charged defects are unfortunately subjected to the spurious 
image interaction, and no supercell size in practice would be sufficient to render this long- 
ranged electrostatic interaction negligible. Various types of corrections have been proposed 
to remove the interactions between the charged defect, its image, and the background charge 
j^j-jd]. On the other hand, in the mixed-basis approach, it is natural to drop the logarith¬ 
mically divergent term of Eq. (|26|) that has similar effect to the cancellation between the 
electrostatic energies from the defects and the uniform compensating background assumed 
in the supercell model. No further corrections are needed in our scheme since only one 
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single isolated charged system rather than an array of the replicated ones is under consid¬ 
eration. Hence, we have developed an alternative promising method to study both neutral 
and charged ID systems with no complications. 


IV. CONCLUSIONS 

In conclusion, we have successfully extended the previous mixed-basis approach to inves¬ 
tigate the electronic structures of one-dimensional systems with plane waves for the periodic 
direction and B-spline sets for the two non-periodic directions. As compared to the existing 
algorithms based upon the conventional supercell model with alternating slab and vacuum 
regions, it is a real space approach along the two non-periodic directions. Therefore, the 
number of the basis functions used to expand the wavefunction is signihcantly reduced 
and the spurious Coulomb interaction between the defect, its images and the compensating 
background charge appeared in the supercell approach can be automatically avoided. 

The new technique has been demonstrated to yield accurate and computationally efficient 
treatments of the inhnite carbon-dimmer chain, graphene nanoribbon, carbon nanotube, and 
the system of the positively-charged chain. It is found that the band structures are all in 
good agreement with those by the popular existing codes, but with a reduced number of basis 
functions. Moreover, no further corrections are needed for the charged case. We have shown 
that the present method is very suitable to investigate either neutral or charged ID materials 
without the need of artificially large supercells or corrections for supercell interactions. 
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Appendix 


1 An 
r (27r)^ 

1 



3*qr 




-d\ 


27r^ 

1 

27r^ 




q\\ + ql " 

g*<?||PCOS0 

q\\ “T dz 


With the identity 


r*27r 


Jn(x) = 


2n 


where Jq is the Bessel function of order 0, 

1 1 


Using another identity 


r TT 


xJo{ax) 

x'^ + k'^ 


dll+Qz 


dx = KQ^ak) [a > 0, Re /c > 0] 


We obtain 


1 _ 1 
r TT 


dq^KoMp) . 


(A.l) 

(A.2) 

(A.3) 

(A.4) 

(A.5) 

(A.6) 

(A.7) 


[1] A. J. Lee, T.-L. Chan, and J. R. Chelikowsky, Phys. Rev. B 89, 075419 (2014). 

[2] S. E. Taylor and F. Bruneval, Phys. Rev. B 84, 075155 (2011), and references therein. 

[3] G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995). 

[4] C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Phys. Rev. B 73, 205119 
(2006). 

[5] G.-W. Li and Y.-G. Ghang, Phys. Rev. B 48, 12032 (1993). 

[6] G.-W. Li and Y.-G. Ghang, Phys. Rev. B 50, 8675 (1994). 

[7] Y.-G. Ghang and G.-W. Li, Gomp. Phys. Gomm. 95, 158 (1996). 

[8] G. Y. Ren, G. S. Hsue and Y.-G. Ghang, Gomp. Phys. Gomm. 188, 94 (2015). 

[9] Garl deBoor, A practical Guide to Splines, (Springer, New York, 1987). 

[10] W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys. Rev. A 37, 307 (1988). 


14 



[11] H. T. Jeng, and C. S. Hsue, Phys. Rev. B 62, 9876 (2000). 

[12] C. Y. Ren, H. T. Jeng, and C. S. Hsue, Phys. Rev. B 66, 125105 (2002). 

[13] D. Vanderbilt, Phys. Rev. B 41, 7982 (1990). 

[14] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 

[15] G. Kresse and J. Furthmiiller, Gomput. Mater. Sci. 6, 15 (1996). 

[16] D. M. Geperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 

[17] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 

[18] L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982). 

[19] K. Laasonen, A. Pasquarello, R. Gar, G. Lee, and D. Vanderbilt, Phys. Rev. B 47, 10142 
(1993). 

[20] http://www.physics.rutgers.edu/~dhv/uspp/. 


15 



FIGURE CAPTIONS 


Fig. 1: Schematic plot of the grid unit containing the singularity at p = p' for the Coulomb 
potential. See text for details. 


Fig. 2: (Color online) (top) Atomic structure of the inhnite carbon-dimmer chain, 
(a) and (b) are the corresponding band structures obtained by the present work and VASP. 


Fig. 3: (Color online) (left) Atomic structure of the armchair graphene nanoribbon, 
a), (b) and (c) are the band structures near Fermi level by the present work, previous work 
8j and VASP, respectively. 


Fig. 4: (Color online) (top) Atomic structure of the (4,0) zigzag carbon nanotube, 
(a) and (b) are the band structures obtained by the present work and VASP. 

Fig. 5: (Color online) The band structures of the positively charged carbon-dimmer 
chain obtained by (a) the present work and (b) VASP with a uniform background charge. 
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